﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Solver
{
  public class Solver
  {
    #region Переменные
    // Число разбиений по x
    int n;
    // Число разбиений по y
    int m;

    //Шаги по пространству
    double h, dt;

    // Число шагов которое совершил метод и максимальное число шагов
    int Smax;
    int S;

    // Достигнутая точность и критерий по точности
    double Eps, Epsmax;

    // Причина останова
    // 1 - по точности
    // 2 - числу шагов
    int clause;

    // вектор-решение,
    // упорядочен слево направо, снизу вверх
    double[] v;

    // Границы области
    double lborderX, lborderY, rborderX, rborderY;

    #endregion

    #region внутренние методы
    void SeidelMethod(int type)
    {
      S = 0;

      double eps_cur = 0;
      double x_new;
      clause = 0;
      Eps = 0;
      while (clause == 0)
      {
        Eps = 0;
        for (int j = 1; j < m; j++)
        {
          for(int i = 1; i < n; i++)
          {
            if ((n / 4 <= i) && (i <= 3 * n / 4) && (m / 4 <= j) && (j <= 3 * m / 4))
            {
              continue;
            }
            // left/right/up/down/top indexes
            int lidx = getIdx(i - 1, j, n);
            int ridx = getIdx(i + 1, j, n);
            int idx = getIdx(i, j, n);
            int uidx = getIdx(i, j + 1, n);
            int didx = getIdx(i, j - 1, n);

            // дополнительный коэффициент
            double A = 2 * ((1.0 / h) * (1.0 / h) + (1.0 / dt) * (1.0 / dt));
            x_new = getF(i,j,type) + (v[lidx] + v[ridx]) / (h * h) 
              + (v[uidx] + v[didx]) / (dt * dt) ;
            x_new = x_new / A;
            if (x_new > 1e3)
            {
              x_new = 0;
            }
            eps_cur = Math.Abs(v[idx] - x_new);
            if (eps_cur > Eps)
            {
              Eps = eps_cur;
            }
            v[idx] = x_new;
          }
        }
        S = S + 1;
        if (Eps < Epsmax)
          clause = 1;
        if (S > Smax)
          clause = 2;
      }
    }
    // Граничные условия основной задачи
    double get0Y(double y)
    {
      return 1 - y * y;
    }
    double getNY(double y)
    {
      return (1 - y * y) * Math.Exp(y);
    }
    double getX0(double x)
    {
      return 1 - x * x;
    }
    double getXM(double x)
    {
      return 1 - x * x;
    }

    // Возмущение (функция f) в задаче Дирихле для уравнения Пуассона
    double getF(int i, int j,int type)
    {
      double answer = -1.0;
      if (type == -1)
      {
        answer = getFtest2(i, j);
      }
      if (type == 0)
      {
        answer = getFtest1(i, j);
      }
      if (type == 1)
      {
        answer = getFMajorTask(i, j);
      }
      return answer;
    }
    double getFMajorTask(int i, int j)
    {
      double x = lborderX + i * h;
      double y = lborderY + j * dt;
      return Math.Abs(x * x - y * y);
    }
    // Возмущение в тестовой задаче 2 типа ( type = -1)
    double getFtest2(int i, int j)
    {
      return 4;
    }
    // Возмущение в тестовой задаче 1 типа ( type = 0)
    double getFtest1(int i , int j)
    {
      double x = lborderX + i * h;
      double y = lborderY + j * dt;
      double arg = 1 - x * x - y * y;
      //4*x^2*e^-x^2-y^2+1-4*e^-x^2-y^2+1+4*y^2*e^-x^2-y^2+1
      return -(4 * x * x * Math.Exp(arg) - 4 * Math.Exp(arg) 
        + 4 * y * y * Math.Exp(arg));
    }

    // Граничное условие тестовой задачи
    double getTestV1(int i, int j)
    {
      double x = lborderX + i * h;
      double y = lborderY + j * dt;
      return Math.Exp(1 - x * x - y * y);
    }
    double getTestV2(int i, int j)
    {
      double x = lborderX + i * h;
      double y = lborderY + j * dt;
      return (1 - x * x - y * y);
    }

    // Вернуть одномерный индекс по двумерному
    public int getIdx(int i, int j, int n)
    {
      return (j * (n + 1) + i);
    }
    #endregion

    #region Внешние методы
    // Последний аргумент тип задачи: -1,0 - тестовые, 1 - основная
    public void prepareParams(int N, int M, int SMAX, double EPSMAX,
      double lx, double ly, double rx, double ry, int type)
    {
      n = N;
      m = M;
      Smax = SMAX;
      Epsmax = EPSMAX;
      lborderX = lx;
      lborderY = ly;
      rborderX = rx;
      rborderY = ry;

      h = (rborderX - lborderX) / n;
      dt = (rborderY - lborderY) / m;

      v = new double[(n + 1) * (m + 1)];
      #region Основная задача
      if (type == 1)
      {
        // Нижняя граница
        for (int i = 0; i <= n; i++)
        {
          int j = 0;
          int idx = getIdx(i, j, n);
          v[idx] = getX0((double)i);
        }

        // Верхняя граница
        for (int i = 0; i <= n; i++)
        {
          int j = m;
          int idx = getIdx(i, j, n);
          v[idx] = getXM((double)i);
        }

        // Левая граница
        for (int j = 0; j <= m; j++)
        {
          int i = 0;
          int idx = getIdx(i, j, n);
          v[idx] = get0Y((double)j);
        }

        // Правая граница
        for (int j = 0; j <= m; j++)
        {
          int i = n;
          int idx = getIdx(i, j, n);
          v[idx] = getNY((double)j);
        }
      }
#endregion
      #region Тестовая задача 1 типа (не полином в решении)
      if (type == 0)
      {

        // Нижняя граница
        for (int i = 0; i <= n; i++)
        {
          int j = 0;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV1(i, j);
        }

        // Верхняя граница
        for (int i = 0; i <= n; i++)
        {
          int j = m;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV1(i, j);
        }

        // Левая граница
        for (int j = 0; j <= m; j++)
        {
          int i = 0;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV1(i, j);
        }

        // Правая граница
        for (int j = 0; j <= m; j++)
        {
          int i = n;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV1(i, j);
        }

        // Внутренние границы

        // Нижняя граница
        for (int i = n / 4; i <= (3 * n / 4); i++)
        {
          int j = m / 4;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV1(i, j);
        }

        // Верхняя граница
        for (int i = n / 4; i <= (3 * n / 4); i++)
        {
          int j = 3 * m / 4;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV1(i, j);
        }

        // Левая граница
        for (int j = m / 4; j <= 3 * m / 4; j++)
        {
          int i = n / 4;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV1(i, j);
        }

        // Правая граница
        for (int j = m / 4; j <= 3 * m / 4; j++)
        {
          int i = 3 * n / 4;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV1(i, j);
        }
      }
#endregion
      #region Тестовая задача 2 типа (полином в решении)
      if (type == -1)
      {

        // Нижняя граница
        for (int i = 0; i <= n; i++)
        {
          int j = 0;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV2(i, j);
        }

        // Верхняя граница
        for (int i = 0; i <= n; i++)
        {
          int j = m;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV2(i, j);
        }

        // Левая граница
        for (int j = 0; j <= m; j++)
        {
          int i = 0;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV2(i, j);
        }

        // Правая граница
        for (int j = 0; j <= m; j++)
        {
          int i = n;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV2(i, j);
        }

        // Внутренние границы

        // Нижняя граница
        for (int i = n / 4; i <= (3 * n / 4); i++)
        {
          int j = m / 4;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV2(i, j);
        }

        // Верхняя граница
        for (int i = n / 4; i <= (3 * n / 4); i++)
        {
          int j = 3 * m / 4;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV2(i, j);
        }

        // Левая граница
        for (int j = m / 4; j <=  3 * m / 4; j++)
        {
          int i = n / 4;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV2(i, j);
        }

        // Правая граница
        for (int j = m / 4; j <= 3 * m / 4; j++)
        {
          int i = 3 * n / 4;
          int idx = getIdx(i, j, n);
          v[idx] = getTestV2(i, j);
        }
      }
      #endregion
    }

    public void solve(int type)
    {
      SeidelMethod(type);
    }

    // узнать причину выхода
    public int getClause()
    {
      return clause;
    }

    // Получить (i,j) компоненту решения
    // 0 <= i <= n
    // 0 <= j <= m
    public double getValueV(int i, int j)
    {
      return v[getIdx(i, j, n)];
    }
    public double getTestSolution(double x, double y, int type)
    {
      double answer = -1.0;
      if (type == 0)
      {
        answer = Math.Exp(1 - x * x - y * y);
      }
      if (type == -1)
      {
        answer = 1 - x * x - y * y;
      }
      return answer;
    }

    // Получить число истраченных шагов
    public int getS()
    {
      return S;
    }

    // Получить достигнутую точность
    public double getEps()
    {
      return Eps;
    }

    #endregion


    // прототипы
    // провести расчёт методом (якоби...)
    // void solve(...);
    // Возвращает текущее решение
    // double returnValue(int i, int j);

    // получить границы
    // 0 <= i <= N
    // Горизонтальные границы
    // getValueX0(double i);
    // getValueXM(double i);
    // 0 <= j <= M
    // Вертикальные границы
    // getValue0Y(double j);
    // getValueNY (double j);
    // Получить значения внутри области
    // getValueF(double x, double y);

  }
}
